01 Finding Concurrent ECOSTRESS and EMIT Data

If running this notebook locally, you will find instructions to set up a compatible environment in the setup folder. If running on the Openscapes 2i2c Cloud Instance for a Workshop, no additional setup is required.

Summary

Both the ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) and the Earth surface Mineral dust source InvesTigation (EMIT) instruments are located on the International Space Station (ISS). Their overlapping fields of view provide an unprecedented opportunity to demonstrate the compounded benefits of working with both datasets. In this notebook we will show how to utilize the earthaccess Python library to find concurrent ECOSTRESS and EMIT data.

Background

The ECOSTRESS insturment is a multispectral thermal imaging radiometer designed to answer three overarching science questions:

• How is the terrestrial biosphere responding to changes in water availability?
• How do changes in diurnal vegetation water stress the global carbon cycle?
• Can agricultural vulnerability be reduced through advanced monitoring of agricultural water consumptive use and improved drought estimation?

The ECOSTRESS mission is answering these questions by accurately measuring the temperature of plants. Plants regulate their temperature by releasing water through tiny pores on their leaves called stomata. If they have sufficient water they can maintain their temperature, but if there is insufficient water, their temperatures rise and this temperature rise can be measured with ECOSTRESS. The images acquired by ECOSTRESS are the most detailed temperature images of the surface ever acquired from space and can be used to measure the temperature of an individual farmers field.

More details about ECOSTRESS and its associated products can be found on the ECOSTRESS website and ECOSTRESS product pages hosted by the Land Processes Distributed Active Archive Center (LP DAAC).

The EMIT instrument is an imaging spectrometer that measures light in visible and infrared wavelengths. These measurements display unique spectral signatures that correspond to the composition on the Earth’s surface. The EMIT mission focuses specifically on mapping the composition of minerals to better understand the effects of mineral dust throughout the Earth system and human populations now and in the future. In addition, the EMIT instrument can be used in other applications, such as mapping of greenhouse gases, snow properties, and water resources.

More details about EMIT and its associated products can be found on the EMIT website and EMIT product pages hosted by the LP DAAC.

Requirements - NASA Earthdata Account - No Python setup requirements if connected to the workshop cloud instance! - Set up Python Environment - See setup_instructions.md in the /setup/ folder

Learning Objectives
- How to use earthaccess to find concurrent EMIT and ECOSTRESS data. - How to export a list of files and download them programatically.

Tutorial Outline

1. Setup
   
2. Searching for Data
   
3. Organizing and Filtering Results
4. Visualizing Intersecting Coverage
5. Creating a list of Asset URLs and Downloading

1. Setup

Import the required Python libraries.

# Import required libraries
import os
import folium
import earthaccess
import warnings
import folium.plugins
import pandas as pd
import geopandas as gpd
import math

from branca.element import Figure
from IPython.display import display
from shapely import geometry
from skimage import io
from datetime import timedelta
from shapely.geometry.polygon import orient
from matplotlib import pyplot as plt

1.2 NASA Earthdata Login Credentials

To download or stream NASA data you will need an Earthdata account, you can create one here. We will use the login function from the earthaccess library for authentication before downloading at the end of the notebook. This function can also be used to create a local .netrc file if it doesn’t exist, or add your login info to an existing .netrc file. If no Earthdata Login credentials are found in the .netrc you’ll be prompted for them. This step is not necessary to conduct searches, but is needed to download or stream data.

2. Search for ECOSTRESS and EMIT Data

Both EMIT and ECOSTRESS products are hosted by the Land Processes Distributed Active Archive Center (LP DAAC). In this example we will use the cloud-hosted EMIT_L2A_RFL and ECOSTRESS_L2T_LSTE products available from the LP DAAC to find data. Any results we find for these products, should be available for other products within the EMIT and ECOSTRESS collections.

To find data we will use the earthaccess Python library. earthaccess searches NASA’s Common Metadata Repository (CMR), a metadata system that catalogs Earth Science data and associated metadata records. The results can then be used to download granules or generate lists granule search result URLs.

Using earthaccess we can search based on the attributes of a granule, which can be thought of as a spatiotemporal scene from an instrument containing multiple assets (eg. Reflectance, Reflectance Uncertainty, Masks for the EMIT L2A Reflectance Collection). We can search using attributes such as collection, acquisition time, and spatial footprint. This process can also be used with other EMIT or ECOSTRESS products, other collections, or different data providers, as well as across multiple catalogs with some modification.

2.1 Define Spatial Region of Interest

For this example, our spatial region of interest (ROI) will be the Carpenteria Salt Marsh. You can learn more about it here: https://ucnrs.org/reserves/carpinteria-salt-marsh-reserve/. If you want to create a geojson polygon for your own ROI, you can do so using this website: https://geojson.io/#map=2/20/0, or you can convert a shapefile to a geojson using some code in the Appendices.

In this example, we elect to search using a polygon rather than a standard bounding box because bounding boxes will have a larger spatial extent, capturing a lot of area we may not be interested in. This becomes more important for searches with larger ROIs than our example here. To search for intersections with a polygon using earthaccess, we need to format our ROI as a counter-clockwise list of coordinate pairs.

Open the geojson file containing a landcover classification of Carpenteria Salt Marsh as a geodataframe, and check the coordinate reference system (CRS) of the data.

polygon = gpd.read_file('../data/landcover.geojson')
polygon.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

The CRS is EPSG:4326 (WGS84), which is also the CRS we want the data in to submit for our search.

Next, lets examine our polygon a bit closer.

polygon

	type	geometry
0	channel	MULTIPOLYGON (((-119.54125 34.40462, -119.5412...
1	salt flat	MULTIPOLYGON (((-119.52907 34.39633, -119.5290...
2	upland	MULTIPOLYGON (((-119.54524 34.40555, -119.5452...
3	pan	MULTIPOLYGON (((-119.52924 34.39675, -119.5292...
4	marsh	MULTIPOLYGON (((-119.54162 34.40421, -119.5416...

We can see this geodataframe consists of multiple classes, each containing a multipolygon within our study site. We need to create an exterior boundary polygon containing these, and make sure the vertices are in counter-clockwise order to submit them in our query. To do this, create a polygon consisting of all the geometries, then calculate the convex hull of the union. This will give us a simple exterior polygon around our full ROI. After that, use the orient function to place our coordinate pairs in counter-clockwise order.

# Merge all Polygon geometries and create external boundary
roi_poly = polygon.unary_union.convex_hull
# Re-order vertices to counter-clockwise
roi_poly = orient(roi_poly, sign=1.0)

We can go ahead and visualize our region of interest and the original landcover polygon. First add a function to help reformat bound box coordinates to work with leaflet notation.

# Function to convert a bounding box for use in leaflet notation

def convert_bounds(bbox, invert_y=False):
    """
    Helper method for changing bounding box representation to leaflet notation

    ``(lon1, lat1, lon2, lat2) -> ((lat1, lon1), (lat2, lon2))``
    """
    x1, y1, x2, y2 = bbox
    if invert_y:
        y1, y2 = y2, y1
    return ((y1, x1), (y2, x2))

Then create a figure using folium.

fig = Figure(width="800px", height="400px")
map1 = folium.Map(tiles='https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}', attr='Google')
fig.add_child(map1)

# Add Convex Hull Polygon
folium.GeoJson(roi_poly,
                name='convex hull',
                ).add_to(map1)

# Add landcover classification geodataframe
polygon.explore(
    "type",
    popup=True,
    categorical=True,
    cmap='Set3',
    style_kwds=dict(opacity=0.7, fillOpacity=0.4),
    name="Carpenteria Salt Marsh Landcover",
    m=map1
)

map1.add_child(folium.LayerControl())
map1.fit_bounds(bounds=convert_bounds(polygon.unary_union.bounds))
display(fig)

Above we can see our region of interest (ROI) and the landcover classification polygon that we opened. We can hover over different areas to see the land cover class.

Lastly we need to convert our polygon to a list of coordinate pairs.

# Set ROI as list of exterior polygon vertices as coordinate pairs
roi = list(roi_poly.exterior.coords)

2.2 Define Collections of Interest

We need to specify which products we want to search for using their short-names. As mentioned above, we will conduct our search using the EMIT Level 2A Reflectance (EMITL2ARFL) and ECOSTRESS Level 2 Tiled Land Surface Temperature and Emmissivity (ECO_L2T_LSTE).

Note: Here we use the Tiled ECOSTRESS LSTE Product. This will also work with the gridded LSTE and the swath; however, the swath product does not have a browse image for the visualization in section 4, and will require additional processing for subsequent analysis.

# Data Collections for our search
collections = ['EMITL2ARFL', 'ECO_L2T_LSTE']

2.3 Define Date Range

For our date range, we’ll look at data collected across the month of April 2023. The date_range can be specified as a pair of dates, start and end (up to, not including).

# Define Date Range
date_range = ('2023-01-01','2023-09-01')

2.4 Searching

Submit a query using earthaccess.

results = earthaccess.search_data(
    short_name=collections,
    polygon=roi,
    temporal=date_range,
    count=500
)

Granules found: 119

len(results)

119

3. Organizing and Filtering Results

As we can see from above, the results object contains a list of objects with metadata and links. We can convert this to a more readable format, a dataframe. In addition, we can make it a geodataframe by taking the spatial metadata and creating a shapely polygon representing the spatial coverage, and further customize which information we want to use from other metadata fields.

First, we define some functions to help us create a shapely object for our geodataframe, and retrieve the specific browse image URLs that we want. By default the browse image selected by earthaccess is the first one in the list, but the ECO_L2_LSTE has several browse images and we want to make sure we retrieve the png file, which is a preview of the LSTE.

# Function to create shapely polygon of spatial coverage
def get_shapely_object(result:earthaccess.results.DataGranule):
    # Get Geometry Keys
    geo = result['umm']['SpatialExtent']['HorizontalSpatialDomain']['Geometry']
    keys = geo.keys()

    if 'BoundingRectangles' in keys:
        bounding_rectangle = geo['BoundingRectangles'][0]
        # Create bbox tuple
        bbox_coords = (bounding_rectangle['WestBoundingCoordinate'],bounding_rectangle['SouthBoundingCoordinate'],
                    bounding_rectangle['EastBoundingCoordinate'],bounding_rectangle['NorthBoundingCoordinate'])
        # Create shapely geometry from bbox
        shape = geometry.box(*bbox_coords, ccw=True)
    elif 'GPolygons' in keys:
        points = geo['GPolygons'][0]['Boundary']['Points']
        # Create shapely geometry from polygons
        shape = geometry.Polygon([[p['Longitude'],p['Latitude']] for p in points])
    else:
         raise ValueError('Provided result does not contain bounding boxes/polygons or is incompatible.')
    return(shape)

# Retrieve png browse image if it exists or first jpg in list of urls
def get_png(result:earthaccess.results.DataGranule):
    https_links = [link for link in result.dataviz_links() if 'https' in link]
    if len(https_links) == 1:
        browse = https_links[0]
    elif len(https_links) == 0:
        browse = 'no browse image'
        warnings.warn(f"There is no browse imagery for {result['umm']['GranuleUR']}.")
    else:
        browse = [png for png in https_links if '.png' in png][0]
    return(browse)

Now that we have our functions we can create a dataframe, then calculate and add our shapely geometries to make a geodataframe. After that, add a column for our browse image urls and print the number of granules in our results, so we can monitor the quantity we are working with a we winnow down to the data we want.

# Create Dataframe of Results Metadata
results_df = pd.json_normalize(results)
# Create shapely polygons for result
geometries = [get_shapely_object(results[index]) for index in results_df.index.to_list()]
# Convert to GeoDataframe
gdf = gpd.GeoDataFrame(results_df, geometry=geometries, crs="EPSG:4326")
# Remove results df, no longer needed
del results_df
# Add browse imagery links
gdf['browse'] = [get_png(granule) for granule in results]
gdf['shortname'] = [result['umm']['CollectionReference']['ShortName'] for result in results]
# Preview GeoDataframe
print(f'{gdf.shape[0]} granules total')

119 granules total

Preview our geodataframe to get an idea what it looks like.

gdf.head()

	size	meta.concept-type	meta.concept-id	meta.revision-id	meta.native-id	meta.provider-id	meta.format	meta.revision-date	umm.TemporalExtent.RangeDateTime.BeginningDateTime	umm.TemporalExtent.RangeDateTime.EndingDateTime	...	umm.Platforms	umm.MetadataSpecification.URL	umm.MetadataSpecification.Name	umm.MetadataSpecification.Version	umm.SpatialExtent.HorizontalSpatialDomain.Geometry.GPolygons	umm.PGEVersionClass.PGEName	umm.CloudCover	geometry	browse	shortname
0	2.49626	granule	G2581836170-LPCLOUD	1	ECOv002_L2T_LSTE_25460_016_11SKU_20230101T1552...	LPCLOUD	application/echo10+xml	2023-01-07T13:31:07.584Z	2023-01-01T15:52:48.650Z	2023-01-01T15:53:40.620Z	...	[{'ShortName': 'ISS', 'Instruments': [{'ShortN...	https://cdn.earthdata.nasa.gov/umm/granule/v1.6.5	UMM-G	1.6.5	NaN	NaN	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
1	8.49543	granule	G2586136993-LPCLOUD	1	ECOv002_L2T_LSTE_25486_006_11SKU_20230103T0744...	LPCLOUD	application/echo10+xml	2023-01-10T22:21:57.631Z	2023-01-03T07:44:22.480Z	2023-01-03T07:45:14.450Z	...	[{'ShortName': 'ISS', 'Instruments': [{'ShortN...	https://cdn.earthdata.nasa.gov/umm/granule/v1.6.5	UMM-G	1.6.5	NaN	NaN	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
2	1.16068	granule	G2591892077-LPCLOUD	1	ECOv002_L2T_LSTE_25547_005_11SKU_20230107T0607...	LPCLOUD	application/echo10+xml	2023-01-18T19:13:09.017Z	2023-01-07T06:07:21.560Z	2023-01-07T06:08:13.530Z	...	[{'ShortName': 'ISS', 'Instruments': [{'ShortN...	https://cdn.earthdata.nasa.gov/umm/granule/v1.6.5	UMM-G	1.6.5	NaN	NaN	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
3	3.42404	granule	G2592560828-LPCLOUD	1	ECOv002_L2T_LSTE_25608_003_11SKU_20230111T0430...	LPCLOUD	application/echo10+xml	2023-01-19T11:23:08.273Z	2023-01-11T04:30:26.510Z	2023-01-11T04:31:18.480Z	...	[{'ShortName': 'ISS', 'Instruments': [{'ShortN...	https://cdn.earthdata.nasa.gov/umm/granule/v1.6.5	UMM-G	1.6.5	NaN	NaN	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
4	14.28010	granule	G2593132604-LPCLOUD	1	ECOv002_L2T_LSTE_25628_015_11SKU_20230112T1150...	LPCLOUD	application/echo10+xml	2023-01-20T05:37:00.651Z	2023-01-12T11:50:28.830Z	2023-01-12T11:51:20.800Z	...	[{'ShortName': 'ISS', 'Instruments': [{'ShortN...	https://cdn.earthdata.nasa.gov/umm/granule/v1.6.5	UMM-G	1.6.5	NaN	NaN	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE

5 rows × 34 columns

There are a lot of columns with data that is not relevant to our goal, so we can drop those. To do that, list the names of colums.

# List Column Names
gdf.columns

Index(['size', 'meta.concept-type', 'meta.concept-id', 'meta.revision-id',
       'meta.native-id', 'meta.provider-id', 'meta.format',
       'meta.revision-date',
       'umm.TemporalExtent.RangeDateTime.BeginningDateTime',
       'umm.TemporalExtent.RangeDateTime.EndingDateTime',
       'umm.OrbitCalculatedSpatialDomains', 'umm.GranuleUR',
       'umm.AdditionalAttributes', 'umm.MeasuredParameters',
       'umm.SpatialExtent.HorizontalSpatialDomain.Geometry.BoundingRectangles',
       'umm.ProviderDates', 'umm.CollectionReference.ShortName',
       'umm.CollectionReference.Version', 'umm.PGEVersionClass.PGEVersion',
       'umm.RelatedUrls', 'umm.DataGranule.DayNightFlag',
       'umm.DataGranule.Identifiers', 'umm.DataGranule.ProductionDateTime',
       'umm.DataGranule.ArchiveAndDistributionInformation', 'umm.Platforms',
       'umm.MetadataSpecification.URL', 'umm.MetadataSpecification.Name',
       'umm.MetadataSpecification.Version',
       'umm.SpatialExtent.HorizontalSpatialDomain.Geometry.GPolygons',
       'umm.PGEVersionClass.PGEName', 'umm.CloudCover', 'geometry', 'browse',
       'shortname'],
      dtype='object')

Now create a list of columns to keep and use it to filter the dataframe.

# Create a list of columns to keep
keep_cols = ['meta.concept-id','meta.native-id', 'umm.TemporalExtent.RangeDateTime.BeginningDateTime','umm.TemporalExtent.RangeDateTime.EndingDateTime','umm.CloudCover','umm.DataGranule.DayNightFlag','geometry','browse', 'shortname']
# Remove unneeded columns
gdf = gdf[gdf.columns.intersection(keep_cols)]
gdf.head()

	meta.concept-id	meta.native-id	umm.TemporalExtent.RangeDateTime.BeginningDateTime	umm.TemporalExtent.RangeDateTime.EndingDateTime	umm.DataGranule.DayNightFlag	umm.CloudCover	geometry	browse	shortname
0	G2581836170-LPCLOUD	ECOv002_L2T_LSTE_25460_016_11SKU_20230101T1552...	2023-01-01T15:52:48.650Z	2023-01-01T15:53:40.620Z	Day	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
1	G2586136993-LPCLOUD	ECOv002_L2T_LSTE_25486_006_11SKU_20230103T0744...	2023-01-03T07:44:22.480Z	2023-01-03T07:45:14.450Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
2	G2591892077-LPCLOUD	ECOv002_L2T_LSTE_25547_005_11SKU_20230107T0607...	2023-01-07T06:07:21.560Z	2023-01-07T06:08:13.530Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
3	G2592560828-LPCLOUD	ECOv002_L2T_LSTE_25608_003_11SKU_20230111T0430...	2023-01-11T04:30:26.510Z	2023-01-11T04:31:18.480Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
4	G2593132604-LPCLOUD	ECOv002_L2T_LSTE_25628_015_11SKU_20230112T1150...	2023-01-12T11:50:28.830Z	2023-01-12T11:51:20.800Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE

This is looking better, but we can make it more readable by renaming our columns.

# Rename some Columns
gdf.rename(columns = {'meta.concept-id':'concept_id','meta.native-id':'granule',
                       'umm.TemporalExtent.RangeDateTime.BeginningDateTime':'start_datetime',
                      'umm.TemporalExtent.RangeDateTime.EndingDateTime':'end_datetime',
                      'umm.CloudCover':'cloud_cover',
                      'umm.DataGranule.DayNightFlag':'day_night'}, inplace=True)
gdf.head()

	concept_id	granule	start_datetime	end_datetime	day_night	cloud_cover	geometry	browse	shortname
0	G2581836170-LPCLOUD	ECOv002_L2T_LSTE_25460_016_11SKU_20230101T1552...	2023-01-01T15:52:48.650Z	2023-01-01T15:53:40.620Z	Day	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
1	G2586136993-LPCLOUD	ECOv002_L2T_LSTE_25486_006_11SKU_20230103T0744...	2023-01-03T07:44:22.480Z	2023-01-03T07:45:14.450Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
2	G2591892077-LPCLOUD	ECOv002_L2T_LSTE_25547_005_11SKU_20230107T0607...	2023-01-07T06:07:21.560Z	2023-01-07T06:08:13.530Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
3	G2592560828-LPCLOUD	ECOv002_L2T_LSTE_25608_003_11SKU_20230111T0430...	2023-01-11T04:30:26.510Z	2023-01-11T04:31:18.480Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE
4	G2593132604-LPCLOUD	ECOv002_L2T_LSTE_25628_015_11SKU_20230112T1150...	2023-01-12T11:50:28.830Z	2023-01-12T11:51:20.800Z	Night	NaN	POLYGON ((-119.06582 34.21003, -119.06582 35.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	ECO_L2T_LSTE

Note: If querying on-premises (not cloud) LP DAAC datasets, the meta.concept-id will not show as xxxxxx-LPCLOUD. For these datasets, the granule name can be retrieved from the umm.DataGranule.Identifiers column.

We can filter using the day/night flag as well, but this step will be unnecessary as we check to ensure all results from ECOSTRESS fall within an hour of resulting EMIT granules.

# gdf = gdf[gdf['day_night'].str.contains('Day')]

Our first step toward filtering the datasets will be to add a column with a datetime.

You may have noticed that the date format is similar for ECOSTRESS and EMIT, but the ECOSTRESS data has an additional fractional seconds. If using the recommended lpdaac_vitals Windows environment, you will need to pass the format='ISO8601'argument to the to_datetime function, like in the commented out line.

# gdf['datetime_obj'] = pd.to_datetime(gdf['start_datetime'])
gdf['datetime_obj'] = pd.to_datetime(gdf['start_datetime'], format='ISO8601')

We can roughly visualize the quantity of results by month at our location using a histogram with 8 bins (Jan up to Sept).

gdf.hist(column='datetime_obj', by='shortname', bins=9, color='green', edgecolor='black', linewidth=1, sharey=True)

array([<Axes: title={'center': 'ECO_L2T_LSTE'}>,
       <Axes: title={'center': 'EMITL2ARFL'}>], dtype=object)

Now we will separate the results into two dataframes, one for ECOTRESS and one for EMIT and print the number of results for each so we can monitor how many granules we’re filtering.

# Suppress Setting with Copy Warning - not applicable in this use case
pd.options.mode.chained_assignment = None  # default='warn'

# Split into two dataframes - ECO and EMIT
eco_gdf = gdf[gdf['granule'].str.contains('ECO')]
emit_gdf = gdf[gdf['granule'].str.contains('EMIT')]
print(f' ECOSTRESS Granules: {eco_gdf.shape[0]} \n EMIT Granules: {emit_gdf.shape[0]}')

 ECOSTRESS Granules: 111 
 EMIT Granules: 8

emit_gdf.head()

	concept_id	granule	start_datetime	end_datetime	day_night	cloud_cover	geometry	browse	shortname	datetime_obj
20	G2631045418-LPCLOUD	EMIT_L2A_RFL_001_20230219T202951_2305013_003	2023-02-19T20:29:51Z	2023-02-19T20:30:14Z	Day	67.0	POLYGON ((-120.04838 35.03646, -120.54523 34.4...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	EMITL2ARFL	2023-02-19 20:29:51+00:00
23	G2631458885-LPCLOUD	EMIT_L2A_RFL_001_20230223T185548_2305412_004	2023-02-23T18:55:48Z	2023-02-23T18:56:00Z	Day	96.0	POLYGON ((-119.95147 35.08163, -120.44685 34.4...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	EMITL2ARFL	2023-02-23 18:55:48+00:00
27	G2631818047-LPCLOUD	EMIT_L2A_RFL_001_20230227T172140_2305811_006	2023-02-27T17:21:40Z	2023-02-27T17:21:52Z	Day	100.0	POLYGON ((-119.41012 34.81982, -119.90931 34.1...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	EMITL2ARFL	2023-02-27 17:21:40+00:00
47	G2667368816-LPCLOUD	EMIT_L2A_RFL_001_20230405T190311_2309513_002	2023-04-05T19:03:11Z	2023-04-05T19:03:23Z	Day	71.0	POLYGON ((-119.87804 34.90691, -120.65480 34.2...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	EMITL2ARFL	2023-04-05 19:03:11+00:00
48	G2667369196-LPCLOUD	EMIT_L2A_RFL_001_20230405T190323_2309513_003	2023-04-05T19:03:23Z	2023-04-05T19:03:35Z	Day	8.0	POLYGON ((-119.22086 35.41499, -120.01088 34.8...	https://data.lpdaac.earthdatacloud.nasa.gov/lp...	EMITL2ARFL	2023-04-05 19:03:23+00:00

We still haven’t filtered the locations where EMIT and ECOSTRESS have data at the same spatial location and time-frame. The EMIT acquisition mask has been added to ECOSTRESS, so in most cases if EMIT is collecting data, so will ECOSTRESS, but there are edge cases where this is not true. To do this we’ll use two filters to catch the edge-cases, and provide an example that can be used with other datasets.

First, since EMIT has a smaller swath width, we can can use a unary union of the spatial coverage present in our geodataframe to filter out ecostress granules that do not overlap with it.

# Subset ECOSTRESS Granules in Geodataframe by intersection with EMIT granules
## Create new column based on intersection with union of EMIT polygons.
eco_gdf['intersects'] = eco_gdf.intersects(emit_gdf.unary_union)
## Apply subsetting
eco_gdf = eco_gdf[eco_gdf['intersects'] == True]
print(f' ECOSTRESS Granules: {eco_gdf.shape[0]} \n EMIT Granules: {emit_gdf.shape[0]}')

 ECOSTRESS Granules: 111 
 EMIT Granules: 8

In this instance, our results aren’t narrowed because our region of interest is smaller than a single EMIT scene. If the spatial ROI was very large, this would be much more unlikely.

Additionally, we want to make sure that data in our results are collected at the same time. For EMIT and ECOSTRESS, the EMIT acquisition mask has been added to the ECOSTRESS mask, meaning that if there is an EMIT scene, there should also be an ECOSTRESS scene acquired at the same time. In practice, however, the timestamps on the scenes can vary slightly. In order to capture this slight variability, we need to use a range instead of a single timestamp to capture concurrent data. To do this, we’ll ensure all ECOSTRESS granule start times fall within 10 minutes of any of the EMIT granules in our results, and vice-versa.

Write a function to evaluate whether these datetime objects fall within 10 minutes of one another using the timedelta function.

# Function to Filter timestamps that do not fall within a time_delta of timestamps from the other acquisition time
def concurrent_match(gdf_a:pd.DataFrame, gdf_b:pd.DataFrame, col_name:str, time_delta:timedelta):
    """
    Cross references dataframes containing a datetime object column and keeps rows in 
    each that fall within the provided timedelta of the other. Acceptable time_delta examples:
    
    months=1
    days=1
    hours=1
    minutes=1
    seconds=1 

    """
    # Match Timestamps from Dataframe A with Time-range of entries in Dataframe B
    # Create empty list
    a_list = []
    # Iterate results for product a based on index values
    for _n in gdf_b.index.to_list():
        # Find where product b is within the window of each product a result
        a_matches = (gdf_a[col_name] > gdf_b[col_name][_n]-time_delta) & (gdf_a[col_name] < gdf_b[col_name][_n]+time_delta)
        # Append list with values
        a_list.append(a_matches)
    # Match Timestamps from Dataframe B with Time-range of entries in Dataframe A
    # Create empy list
    b_list =[]
    for _m in gdf_a.index.to_list():
        # Find where product a is within the window of each product b result
        b_matches = (gdf_b[col_name] > gdf_a[col_name][_m]-time_delta) &  (gdf_b[col_name] < gdf_a[col_name][_m]+time_delta)
        # Append list with values
        b_list.append(b_matches)
    # Filter Original Dataframes by summing list of bools, 0 = outside of all time-ranges
    a_filtered = gdf_a.loc[sum(a_list) > 0]
    b_filtered = gdf_b.loc[sum(b_list) > 0]
    return(a_filtered, b_filtered)

Now run our function.

eco_gdf, emit_gdf = concurrent_match(eco_gdf,emit_gdf, col_name='datetime_obj',time_delta=timedelta(minutes=10))
print(f' ECOSTRESS Granules: {eco_gdf.shape[0]} \n EMIT Granules: {emit_gdf.shape[0]}')

 ECOSTRESS Granules: 5 
 EMIT Granules: 6

4. Visualizing Intersecting Coverage

Now that we have geodataframes containing some concurrent data, we can visualize them on a map using folium. It’s often difficult to visualize a large time-series of scenes, so we’ve included an example in Appendix A1 on how to filter to a single day.

# Plot Using Folium

# Create Figure and Select Background Tiles
fig = Figure(width="1100px", height="550px")
map1 = folium.Map(tiles='https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}', attr='Google')
fig.add_child(map1)

# Plot STAC ECOSTRESS Results - note we must drop the datetime_obj columns for this to work
eco_gdf.drop(columns=['datetime_obj']).explore(
    "granule",
    categorical=True,
    tooltip=[
        "granule",
        "start_datetime",
        "cloud_cover",
    ],
    popup=True,
    style_kwds=dict(fillOpacity=0.1, width=2),
    name="ECOSTRESS",
    m=map1,
)

# Plot STAC EMITL2ARFL Results - note we must drop the datetime_obj columns for this to work
emit_gdf.drop(columns=['datetime_obj']).explore(
    "granule",
    categorical=True,
    tooltip=[
        "granule",
        "start_datetime",
        "cloud_cover",
    ],
    popup=True,
    style_kwds=dict(fillOpacity=0.1, width=2),
    name="EMIT",
    m=map1,
)

# ECOSTRESS Browse Images - Comment out to remove
for _n in eco_gdf.index.to_list():
    folium.raster_layers.ImageOverlay(
        image=eco_gdf['browse'][_n],
        name=eco_gdf['granule'][_n],
        bounds=[[eco_gdf.bounds['miny'][_n], eco_gdf.bounds['minx'][_n]], [eco_gdf.bounds['maxy'][_n], eco_gdf.bounds['maxx'][_n]]],
        interactive=False,
        cross_origin=False,
        opacity=0.75,
        zindex=1,
        ).add_to(map1)

# Plot Region of Interest
polygon.explore(
    popup=False,
    style_kwds=dict(fillOpacity=0.1, width=2),
    name="Region of Interest",
    m=map1
)

map1.fit_bounds(bounds=convert_bounds(gdf.unary_union.bounds))
map1.add_child(folium.LayerControl())
display(fig)

In the figure above, you can zoom in and out, click and drag to reposition the legend, and add or remove layers using the layer control in the top right. Notice that since we’re using the tiled ecostress product, we have 2 overlapping tiles at our ROI. You can visualize the tiles by adding or removing the layers.

4.2 Previewing EMIT Browse Imagery

The EMIT browse imagery is not orthorectified, so it can’t be visualized on a plot like the ECOSTRESS browse imagery. To get an idea what scenes look like we can plot them in a grid using matplotlib.

Note: The black space is indicative of onboard cloud masking that occurs before data is downlinked from the ISS.

cols = 3
rows = math.ceil(len(emit_gdf)/cols)
fig, ax = plt.subplots(rows, cols, figsize=(20,20))
ax = ax.flatten()

for _n, index in enumerate(emit_gdf.index.to_list()):
    img = io.imread(emit_gdf['browse'][index])
    ax[_n].imshow(img)
    ax[_n].set_title(f"Index: {index} - {emit_gdf['granule'][index]}")
    ax[_n].axis('off')
plt.tight_layout()
plt.show

<function matplotlib.pyplot.show(close=None, block=None)>

4.3 Further Filtering

We can see that some of these granules likely won’t work because of the large amount of cloud cover, we can use a list of these to filter them out. Make a list of indexes to filter out.

# Bad granule list
bad_granules = [27,74,87]

Filter out the bad granules.

emit_gdf = emit_gdf[~emit_gdf.index.isin(bad_granules)]

Now that we’ve narrowed our EMIT results we can again filter the ecostress granules based on their concurrency with our filtered EMIT granules.

eco_gdf, emit_gdf = concurrent_match(eco_gdf,emit_gdf, col_name='datetime_obj',time_delta=timedelta(hours=1))
print(f' ECOSTRESS Granules: {eco_gdf.shape[0]} \n EMIT Granules: {emit_gdf.shape[0]}')

 ECOSTRESS Granules: 2 
 EMIT Granules: 3

5. Generating a list of URLs and downloading data

Creating a list of results URLs will include all of these assets, so if we only want a subset we need an additional filter to keep the specific assets we want.

If you look back, you can see we kept the same indexing throughout the notebook. This enables us to simply subset the earthaccess results object to retrieve the results we want.

Create a list of index values to keep.

keep_granules = eco_gdf.index.to_list()+emit_gdf.index.to_list()
keep_granules.sort()

Filter the results list.

filtered_results = [result for i, result in enumerate(results) if i in keep_granules]

Now we can download all of the associated assets, or retrieve the URLS and further filter them to specifically what we want.

First, log into Earthdata using the login function from the earthaccess library. The persist=True argument will create a local .netrc file if it doesn’t exist, or add your login info to an existing .netrc file. If no Earthdata Login credentials are found in the .netrc you’ll be prompted for them. As mentioned in section 1.2, this step is not necessary to conduct searches, but is needed to download or stream data.

earthaccess.login(persist=True)

EARTHDATA_USERNAME and EARTHDATA_PASSWORD are not set in the current environment, try setting them or use a different strategy (netrc, interactive)
You're now authenticated with NASA Earthdata Login
Using token with expiration date: 12/24/2023
Using .netrc file for EDL

<earthaccess.auth.Auth at 0x248568a2740>

Now we can download all assets using the following cell.

# # Download All Assets for Granules in Filtered Results
# earthaccess.download(filtered_results, '../data/')

Or we can create a list of URLs and use that to further refine which files we download.

# Retrieve URLS for Assets
results_urls = [granule.data_links() for granule in filtered_results]

Granules often have several assets associated with them, for example, ECO_L2T_LSTE has several assets: - Water Mask (water) - Cloud Mask (cloud) - Quality (QC) - Land Surface Temperature (LST) - Land Surface Temperature Error (LST_err) - Wide Band Emissivity (EmisWB) - Height (height)

The results list we just generated contains URLs to all of these files. We can further filter our results list using string matching to remove unwanted assets.

Create a list of strings and enumerate through our results_url list to filter out unwanted assets.

filtered_asset_links = []
# Pick Desired Assets (leave _ on RFL to distinguish from RFLUNC, LST. to distinguish from LST_err)
desired_assets = ['RFL_', 'LST.'] # Add more or do individually for reflectance, reflectance uncertainty, or mask
# Step through each sublist (granule) and filter based on desired assets.
for n, granule in enumerate(results_urls):
    for url in granule: 
        asset_name = url.split('/')[-1]
        if any(asset in asset_name for asset in desired_assets):
            filtered_asset_links.append(url)
filtered_asset_links

['https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/ECO_L2T_LSTE.002/ECOv002_L2T_LSTE_26223_012_11SKU_20230219T202943_0710_01/ECOv002_L2T_LSTE_26223_012_11SKU_20230219T202943_0710_01_LST.tif',
 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20230219T202951_2305013_003/EMIT_L2A_RFL_001_20230219T202951_2305013_003.nc',
 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/ECO_L2T_LSTE.002/ECOv002_L2T_LSTE_26921_001_11SKU_20230405T190258_0710_01/ECOv002_L2T_LSTE_26921_001_11SKU_20230405T190258_0710_01_LST.tif',
 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20230405T190311_2309513_002/EMIT_L2A_RFL_001_20230405T190311_2309513_002.nc',
 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20230405T190323_2309513_003/EMIT_L2A_RFL_001_20230405T190323_2309513_003.nc']

Uncomment the cell below (select all, then ctrl+/) and download the data that we’ve filtered.

# # Get requests https Session using Earthdata Login Info
# fs = earthaccess.get_requests_https_session()
# # Retrieve granule asset ID from URL (to maintain existing naming convention)
# for url in filtered_asset_links:
#     granule_asset_id = url.split('/')[-1]
#     # Define Local Filepath
#     fp = f'../data/{granule_asset_id}'
#     # Download the Granule Asset if it doesn't exist
#     if not os.path.isfile(fp):
#         with fs.get(url,stream=True) as src:
#             with open(fp,'wb') as dst:
#                 for chunk in src.iter_content(chunk_size=64*1024*1024):
#                     dst.write(chunk)

Congratulations, now you you have downloaded concurrent data from the ECOSTRESS and EMIT instruments on the ISS.

Contact Info:

Email: LPDAAC@usgs.gov
Voice: +1-866-573-3222
Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
Website: https://lpdaac.usgs.gov/
Date last modified: 10-12-2023

¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.

Appendices

These contain some extra code that may be useful when performing a similar workflow.

A1. Further Limiting Search for Visualization Purposes

A large quantity of results may be difficult to understand when mapping with folium. We can create a subset that falls within a single day. First add another column of dates only, then find the unique dates.

# eco_gdf['date'] = eco_gdf['start_datetime'].str.split('T').str[0]
# emit_gdf['date'] = emit_gdf['start_datetime'].str.split('T').str[0]
# emit_gdf['date'].unique()

Filter both sets of results using a single date.

# single_day_eco = eco_gdf#[eco_gdf['date'] == '2023-04-23']
# single_day_emit = emit_gdf#[emit_gdf['date'] == '2023-04-23']
# print(f' ECOSTRESS Granules: {single_day_eco.shape[0]} \n EMIT Granules: {single_day_emit.shape[0]}')

A2. Convert Shapefile to GeoJSON

We can convert a shapefile to a geojson using the following cell. Note that we need to reorder the polygon external vertices so we can submit them as a list of points for our search.

# # Use Sedgwick Reserve Shapefile
# # Open Shapefile
# polygon = gpd.read_file('../data/Sedgwick_Boundary/Sedgwick_Boundary.shp').to_crs("EPSG:4326")
# # Reorder vertices into Counter-clockwise order
# polygon.geometry[0] = orient(polygon.geometry[0], sign=1.0)
# # Save as a geojson (not necessary)
# polygon.to_file('../data/sedgwick_boundary_epsg4326.geojson', driver='GeoJSON')